#log# Automatic Logger file. *** THIS MUST BE THE FIRST LINE ***
#log# DO NOT CHANGE THIS LINE OR THE TWO BELOW
import astro.llsvel
import numpy as np
import pylab as pl

def calc_hist(bins,data):
    """ Given an array of bin edge positions and a data set, returns a
    histogram (number of data points falling in each bin). Note the
    bins must be continuous, and bins is the edge values of the
    bins. i.e. len(bins) = len(data) + 1.  The bins can be different
    widths. If a data value falls on a bin edge, it is placed in the
    lower bin.

    To plot a histogram using the output from this program, use e.g.
    
    >>> binvalues = calc_hist(bins,data)                   # doctest: +SKIP
    >>> pylab.bar(bins[:-1],binvalues,width=0.98*binwidth) # doctest: +SKIP
    """
    counts = []
    for i in xrange(len(bins)-1):
        num = 0.0
        for value in data:
            if bins[i] < value <= bins[i+1]:
                num += 1.0
        counts.append(num)
    return np.array(counts)

def hist_plot(bins,loz,hiz,nsys_loz,nsys_hiz,
              name='hist_med_50kms',title='Median'):
    f = pl.gcf()
    pl.clf()
    # First normalised by number of systems
    f.add_subplot(211)
    pl.bar(bins[:-1],hiz/nsys_hiz,width=binwidth*0.98,alpha=0.5,fc='r')
    pl.bar(bins[:-1],loz/nsys_loz,width=binwidth*0.98,alpha=0.5,fc='b')
    pl.xlabel('Velocity (km/s)')
    pl.ylabel('No. components /\n Total no. of systems')
    pl.title(title)
    # Then by number of components
    f.add_subplot(212)
    pl.bar(bins[:-1],hiz/len(hiz),width=binwidth*0.98,alpha=0.5,fc='r')
    pl.bar(bins[:-1],loz/len(loz),width=binwidth*0.98,alpha=0.5,fc='b')
    pl.xlabel('Velocity (km/s)')
    pl.ylabel('No. components /\nTotal no. of components')
    pl.savefig(name)

hiz_med,hiv_med = astro.llsvel.getcomps('in_lls_vel',centre='median')
hiz_mean,hiv_mean = astro.llsvel.getcomps('in_lls_vel',centre='mean')

loz_med,lov_med = astro.llsvel.getcomps('in_lls_vel_loz',centre='median')
loz_mean,lov_mean = astro.llsvel.getcomps('in_lls_vel_loz',centre='mean')

nloz = len(loz_med)
nhiz = len(hiz_med)

binwidth = 50.
bins = np.arange(0,2000,binwidth) - (2000.-binwidth)/2.

# just use means for now.

countslo = calc_hist(bins,lov_mean)
countshi = calc_hist(bins,hiv_mean)

hist_plot(bins,countslo,countshi,nloz,nhiz,name='hist_mean_50kms',title='Mean')

countslo = calc_hist(bins,lov_med)
countshi = calc_hist(bins,hiv_med)

hist_plot(bins,countslo,countshi,nloz,nhiz,name='hist_med_50kms',title='Median')

binwidth = 100.
bins = np.arange(0,2000,binwidth) - (2000. - binwidth)/2.

countslo = calc_hist(bins,lov_mean)
countshi = calc_hist(bins,hiv_mean)

hist_plot(bins,countslo,countshi,nloz,nhiz,name='hist_mean_100kms',title='Mean')

countslo = calc_hist(bins,lov_med)
countshi = calc_hist(bins,hiv_med)

hist_plot(bins,countslo,countshi,nloz,nhiz,name='hist_med_100kms',title='Median')
